// Statistical catch-at-age for the msur population 
//==+==+==+==+==+==+==+==+==

GLOBALS_SECTION
  #include <logLike.scape>
   
TOP_OF_MAIN_SECTION
  arrmblsize = 50000000;
  gradient_structure::set_NUM_DEPENDENT_VARIABLES(5000);
  gradient_structure::set_MAX_NVAR_OFFSET(1000);
  gradient_structure::set_GRADSTACK_BUFFER_SIZE(400000);
  gradient_structure::set_CMPDIF_BUFFER_SIZE(2000000);

DATA_SECTION
  init_int debug        // flag -- 0
  init_int ncoh		// number of cohorts -- 17
  init_int nages   	// maximum number of month -- 48
  init_int nper         // number of periods in month -- 186
  init_int nyrs         // number of years -- 16
  init_int nlen         // number of length intevals -- 31
  init_number stlen   	// first length interval -- 3
  init_number steplen  	// size step -- 0.5
  
 // Limits lower and upper limits por ragged array abundance matrix  
  init_ivector low(1,ncoh)
  init_ivector up(1,ncoh)
  
 // time period - year relationship
  init_ivector year(1,nper)
  init_ivector meses(1,nper)
  init_imatrix edades(1,ncoh,low,up)

 // ABUNDANCE INDICES
  init_int nobs_Catch		// mensual -- 186
  init_int nobs_CPUE		// mensual -- 168	
  init_int nobs_Reclas		//RECLAS -- 11
  init_int nobs_Pelaces		//PELACES -- 8

 // LENGTH STRUCTURE DATA
  init_int nobs_CL              //CL por mes -- 160 
  init_int nobs_CL_Reclas       //CL Reclas -- 11
  init_int nobs_CL_Pelaces      //CL Pelaces -- 8

 // ABUNDANCE INDICES position
  init_ivector yrs_catch(1,nobs_Catch) //indx meses con capturas -- 186
  init_ivector yrs_CPUE(1,nobs_CPUE) //indx meses con cpue -- 160
  init_ivector yrs_Reclas(1,nobs_Reclas) //indx Reclas meses -- 11 
  init_ivector yrs_Pelaces(1,nobs_Pelaces) //indx Pelaces meses -- 8

// CATCH
  init_vector obs_Catch(1,nobs_Catch)  	 //Catch in tons by month

// ABUNDANCE INDICES
  init_vector obs_CPUE(1,nobs_CPUE) //CPUE 
  init_vector obs_Reclas(1,nobs_Reclas) //Reclas 
  init_vector obs_Pelaces(1,nobs_Pelaces) //Pelaces 
 
 // LENGTH DATA  
  init_ivector yrs_CL(1,nobs_CL) //indx length-structure meses -- 160
  init_ivector yrs_CL_Reclas(1,nobs_CL_Reclas) //indx length-stru -- 11 
  init_ivector yrs_CL_Pelaces(1,nobs_CL_Pelaces) //indx length-stru -- 8 
 
 // BIOLOGICAL DATA
 // GROWTH DATA (BY Cohort) 1:Linf; 2:k; 3:to; 4:A; 5:ts; 6:sd1; 7:sdn
  init_matrix g(1,ncoh,1,7) 

 // MATURITY OGIVE at age(ONE FOR ALL YEARs) 48 observations
  init_number A50_m 
  init_number rango_m
  
 // WEIGHT-AT-LENGTH (BY MONTH)
  init_matrix w_length(1,nper,1,2) // aL^b por periodo -- 186

 // PROPORTIONS LENGTH DATA  
  init_matrix obs_CL(1,nobs_CL,1,nlen) //length-struc fishery -- 160 
  init_matrix obs_CL_Reclas(1,nobs_CL_Reclas,1,nlen) //length-struc reclas -- 11
  init_matrix obs_CL_Pelaces(1,nobs_CL_Pelaces,1,nlen) //length-stru pelaces -- 8 

 //temporal and cohort indexes
  init_ivector ncoh_alive(1,nper);  
  init_imatrix ncoh_which(1,nper,1,ncoh_alive);

 // flags used to determine which likelihoods are being used
  init_int flag_Catch;  
  init_int flag_CPUE;  
  init_int flag_Reclas;  
  init_int flag_Pelaces;  
  init_int flag_CL;  
  init_int flag_CL_Reclas;  
  init_int flag_CL_Pelaces;

 // Lee controles de MOSARDINA
 // --------------------------------------------------------------------------------
 !! ad_comm::change_datafile_name("mosardina_control.dat");

 // CV ABUNDANCE INDICES 
  init_number cvv_Catch  //CV Catch
  init_number cvv_CPUE //CV CPUE 
  init_number cvv_Reclas //CV Reclas 
  init_number cvv_Pelaces //CV Pelaces 

  vector cv_Catch(1,nobs_Catch)  //CV Catch
  !! cv_Catch = cvv_Catch;
  vector cv_CPUE(1,nobs_CPUE) //CV CPUE
  !! cv_CPUE = cvv_CPUE; 
  vector cv_Reclas(1,nobs_Reclas) //CV Reclas
  !! cv_Reclas = cvv_Reclas; 
  vector cv_Pelaces(1,nobs_Pelaces) //CV Pelaces
  !! cv_Pelaces = cvv_Pelaces;

 // SAMPLE SIZE PROPORTIONS LENGTH DATA  
  init_int nn_CL //CPUE 
  init_int nn_CL_Reclas //Reclas 
  init_int nn_CL_Pelaces //Pelaces 

  ivector n_CL(1,nobs_CL) //CPUE
  !! n_CL = nn_CL; 
  ivector n_CL_Reclas(1,nobs_CL_Reclas) //Reclas
  !! n_CL_Reclas = nn_CL_Reclas; 
  ivector n_CL_Pelaces(1,nobs_CL_Pelaces) //Pelaces
  !! n_CL_Pelaces = nn_CL_Pelaces;

 // Bound and phase to all parameters
  init_vector M_dat(1,3)              // Mortality y
  init_vector Afull_dat(1,3)          // A50 fishery
  init_vector varL_dat(1,3)           // sdt^2 left fishery
  init_vector Reclas_A_50_dat(1,3)    // A50 Reclas
  init_vector Reclas_ran_dat(1,3)     // A95-A50 Reclas
  init_vector Pelaces_A_50_dat(1,3)   // A50 Pelaces  
  init_vector Pelaces_ran_dat(1,3)    // A95-A50 Pelaces

  init_vector Ft_dat(1,3) 	      // Anual Fishing Mortality
  init_vector mu1F_dat(1,3) 
  init_vector mu2F_dat(1,3) 
  init_vector sd1F_dat(1,3) 
  init_vector sd2F_dat(1,3) 
  init_vector pF_dat(1,3)             // 
  init_vector log_q_dat(1,3);  	      // catchability fishery
  init_vector Reclas_log_q_dat(1,3);  // proportional observ. Reclas
  init_vector Pelaces_log_q_dat(1,3); // proportional observ. Pelaces
  init_vector log_Rec_dat(1,3);       // Recruitment

 // Vectors and Matrices data-dependents
 // 3d vector for age-length conversion
  3darray Age_length(1,ncoh,1,nages,1,nlen);
  matrix W_length(1,nper,1,nlen);
  matrix mu(1,ncoh,1,nages);
  matrix sd(1,ncoh,1,nages);
  vector length_(1,nlen);
  vector matur(1,nages);
  imatrix mes(1,ncoh,1,nages);
  number pi
  number lb
  number ub
  number suma
  int mes_cal
  int edad
  int phz
  int cont
  int yr
  int id
  int coh
  int LikeType
  int j
  int i
  int k
  int l

PARAMETER_SECTION

 // Natural mortality
 LOCAL_CALCS
   phz=M_dat(1);
   lb=M_dat(2);
   ub=M_dat(3);
 END_CALCS
  init_bounded_number M(lb,ub,phz)

 //Fleet selectivity 
 LOCAL_CALCS
   phz=Afull_dat(1);
   lb=Afull_dat(2);
   ub=Afull_dat(3);
 END_CALCS
  init_bounded_number Afull(lb,ub,phz)  
 
 LOCAL_CALCS
   phz=varL_dat(1);
   lb=varL_dat(2);
   ub=varL_dat(3);
 END_CALCS
  init_bounded_number varL(lb,ub,phz)  
 
 //Reclas selectivity  
 LOCAL_CALCS
   phz=Reclas_A_50_dat(1);
   lb=Reclas_A_50_dat(2);
   ub=Reclas_A_50_dat(3);
 END_CALCS
  init_bounded_number Reclas_A_50(lb,ub,phz)  

 LOCAL_CALCS
   phz=Reclas_ran_dat(1);
   lb=Reclas_ran_dat(2);
   ub=Reclas_ran_dat(3);
 END_CALCS
  init_bounded_number Reclas_ran(lb,ub,phz)  

 //Pelaces selectivity  
 LOCAL_CALCS
   phz=Pelaces_A_50_dat(1);
   lb=Pelaces_A_50_dat(2);
   ub=Pelaces_A_50_dat(3);
 END_CALCS
  init_bounded_number Pelaces_A_50(lb,ub,phz)  
 
 LOCAL_CALCS
   phz=Pelaces_ran_dat(1);
   lb=Pelaces_ran_dat(2);
   ub=Pelaces_ran_dat(3);
 END_CALCS
  init_bounded_number Pelaces_ran(lb,ub,phz)  

 //Fleet catchability 
 LOCAL_CALCS
   phz=log_q_dat(1);
   lb=log_q_dat(2);
   ub=log_q_dat(3);
 END_CALCS
  init_bounded_number log_q(lb,ub,phz)  

 //Reclas catchability  
 LOCAL_CALCS
   phz=Reclas_log_q_dat(1);
   lb=Reclas_log_q_dat(2);
   ub=Reclas_log_q_dat(3);
 END_CALCS
  init_bounded_number Reclas_log_q(lb,ub,phz)  

 //Pelaces catchability   
 LOCAL_CALCS
   phz=Pelaces_log_q_dat(1);
   lb=Pelaces_log_q_dat(2);
   ub=Pelaces_log_q_dat(3);
 END_CALCS
  init_bounded_number Pelaces_log_q(lb,ub,phz)  
 
 // Recruitment
 LOCAL_CALCS
      phz= log_Rec_dat(1);
      lb = log_Rec_dat(2);
      ub = log_Rec_dat(3);
 END_CALCS
  init_bounded_vector log_Rec(1,ncoh,lb,ub,phz)
 
 //Fishing mortality (anual and seasonal distribution - double normal)
 LOCAL_CALCS
      phz= Ft_dat(1);
      lb = Ft_dat(2);
      ub = Ft_dat(3);
 END_CALCS
  init_bounded_vector log_Ft(1,nyrs,lb,ub,phz)
 
 LOCAL_CALCS
      phz= mu1F_dat(1);
      lb = mu1F_dat(2);
      ub = mu1F_dat(3);
 END_CALCS
  init_bounded_vector mu1F(1,nyrs,lb,ub,phz)
 
 LOCAL_CALCS
      phz= mu2F_dat(1);
      lb = mu2F_dat(2);
      ub = mu2F_dat(3);
 END_CALCS
  init_bounded_vector mu2F(1,nyrs,lb,ub,phz)

 LOCAL_CALCS
      phz= sd1F_dat(1);
      lb = sd1F_dat(2);
      ub = sd1F_dat(3);
 END_CALCS
  init_bounded_vector sd1F(1,nyrs,lb,ub,phz)
 
 LOCAL_CALCS
      phz= sd2F_dat(1);
      lb = sd2F_dat(2);
      ub = sd2F_dat(3);
 END_CALCS
  init_bounded_vector sd2F(1,nyrs,lb,ub,phz)
 
 LOCAL_CALCS
      phz= pF_dat(1);
      lb = pF_dat(2);
      ub = pF_dat(3);
 END_CALCS
  init_bounded_vector pF(1,nyrs,lb,ub,phz)
  
    
  //init_number dummy
   
 //*********************   
 // Estimated variables   
 //*********************

 //Catch and CPUE
  vector Catch(1,nper);
  vector pCatch(1,nobs_Catch);
  vector Catch_LL(1,nobs_Catch);
  vector Catch_residuals(1,nobs_Catch);
    
  vector CPUE(1,nper);
  vector pCPUE(1,nobs_CPUE);
  vector CPUE_LL(1,nobs_CPUE);
  vector CPUE_residuals(1,nobs_CPUE);

 //Biomasses
  vector BT(1,nper);         //total biomass
  vector BV(1,nper);         //vulnerable biomass
  vector SB(1,nper);         //spawning biomass
  vector MA(1,nper);         //monthly abundance
  vector RB(1,nper);         //Recruits biomass        
  
 //Abundances and catch at length 
  matrix N(1,ncoh,low,up);
  matrix C(1,ncoh,low,up);
  matrix CL(1,nper,1,nlen); // catch at length
  matrix predCL(1,nobs_CL,1,nlen); 
  matrix CL_LL(1,nobs_CL,1,nlen);
  matrix CL_residuals(1,nobs_CL,1,nlen);
  vector CL_obs_stdev(1,nobs_CL);
  
  matrix tNatL(1,nper,1,nlen); // total abundance at length
  matrix vNatL(1,nper,1,nlen); // spawning abundance at length
  matrix sNatL(1,nper,1,nlen); // vulnerable abundance at length
  matrix vNatL_Reclas_CL(1,nobs_CL_Reclas,1,nlen);  
  matrix vNatL_Reclas_Index(1,nobs_Reclas,1,nlen);  
  matrix vNatL_Pelaces_CL(1,nobs_CL_Pelaces,1,nlen);  
  matrix vNatL_Pelaces_Index(1,nobs_Pelaces,1,nlen);  

  matrix pCL(1,nper,1,nlen); // catch at length

  matrix CL_Reclas(1,nobs_CL_Reclas,1,nlen);
  matrix CL_Reclas_LL(1,nobs_CL_Reclas,1,nlen);
  matrix CL_Reclas_residuals(1,nobs_CL_Reclas,1,nlen);

  matrix CL_Pelaces(1,nobs_CL_Pelaces,1,nlen);
  matrix CL_Pelaces_LL(1,nobs_CL_Pelaces,1,nlen);
  matrix CL_Pelaces_residuals(1,nobs_CL_Pelaces,1,nlen);
  

 //Acoustic indexes
  vector Index_Reclas(1,nobs_Reclas); 
  vector Reclas_LL(1,nobs_Reclas);
  vector Reclas_residuals(1,nobs_Reclas);
  
  
  vector Index_Pelaces(1,nobs_Pelaces); 
  vector Pelaces_LL(1,nobs_Pelaces);
  vector Pelaces_residuals(1,nobs_Pelaces); 
 
 //maturity and selectivities
  vector sel(1,nages);
  vector sel_Reclas(1,nages);
  vector sel_Pelaces(1,nages);
  matrix Fm(1,nyrs,1,12);
  matrix Fm_eff(1,nyrs,1,12);
  number c;
  number c2;
  vector Ft;
  
  sdreport_number f_prof_cl;
  // likeprof_number f_prof_cl;
    
 //init_number dummy(1);

  objective_function_value f


PRELIMINARY_CALCS_SECTION
  pi=3.141592653590;

 //length classes
  length_(1)=stlen;
  for (j=2;j<=nlen;j++)
    length_(j)=length_(j-1)+steplen;
   
 //get_W_length
  for (i=1;i<=nper;i++)
    for (j=1;j<=nlen;j++)
       W_length(i,j)=(w_length(i,1)*pow(length_(j),w_length(i,2)))/1e6; 
       
 
 //get Age-length key 
   for (i=1;i<=ncoh; i++)
      {
        for (j=1;j<=nages; j++)
           {
             mes(i,j)=(i-1)*12+j; //esta matriz puede ser ingresada desde el archivo de datos para incoporar mayor flexibilidad  
             sd(i,j)=g(i,6)-((g(i,6)-g(i,7))/47)*(j-1);
             mu(i,j)= g(i,1)*( 1.0-exp( -g(i,2)* (j/12.0-g(i,3)) -((g(i,4)*g(i,2))/(2*pi)) *(sin( 2*pi* (j/12.0- g(i,5))) - sin(2*pi*(g(i,3)-g(i,5))))));

      	     for (int k=1;k<=nlen; k++)
                {
                  Age_length(i,j,k)= (1/sqrt(2*pi*(square(sd(i,j)))))*exp(-0.5*(square(length_(k)-mu(i,j)))/(square(sd(i,j))));
                  suma=suma+Age_length(i,j,k);
                }
                  Age_length(i,j)=Age_length(i,j)/suma;
             suma=0.0;
           }
      }

      
PROCEDURE_SECTION
  CL.initialize();
  tNatL.initialize();
  vNatL.initialize();
  sNatL.initialize();
  CPUE.initialize();
  Catch.initialize();  
  BT.initialize();
  BV.initialize();
  SB.initialize();
  vNatL_Reclas_CL.initialize();  
  vNatL_Reclas_Index.initialize();  
  vNatL_Pelaces_CL.initialize();  
  vNatL_Pelaces_Index.initialize();  
  Index_Reclas.initialize(); 
  Index_Pelaces.initialize();   
  CL_Reclas.initialize();   
  CL_Pelaces.initialize();   

  f=0; 			// Reset total objective function 
     
  get_selectivities(); 
  get_Fm();
  get_dynamic();
  get_Likelihoods();
  get_DFLs();
  get_residuals();
   
FUNCTION get_selectivities

 //Fleet selectivity 
  for (j=1;j<=nages;j++)
	  {
       if ((Afull)>j)
         {
           sel(j)=exp(-1.0*(square(j-(Afull)))/varL);             
         }
       else
         {
           sel(j)=1.0; 
         }
     }
 
 //Reclas and Pelaces selectivities 
    for (j=1;j<=nages; j++)
	    {
       sel_Reclas(j)=1.0/(1.0+exp(-(log(19.0)*(j-(Reclas_A_50)))/(Reclas_ran)));
       sel_Pelaces(j)=1.0/(1.0+exp(-(log(19.0)*(j-(Pelaces_A_50)))/(Pelaces_ran)));
       matur(j)=1.0/(1.0+exp(-(log(19.0)*(j-(A50_m)))/(rango_m)));
      }


FUNCTION get_Fm
 //Seasonal fishing mortality distributions
    for (i=1;i<=nyrs; i++)
       {
         c=0.0;
         c2=0.0; 
         for (j=1;j<=12; j++)
            {
              Fm(i,j)=pF(i)*(1.0/sqrt(2.0*pi*(square(sd1F(i)))))*exp(-(square(j-mu1F(i)))/(2*square(sd1F(i))))  +  (1.0 - pF(i))*(1.0/sqrt(2.0*pi*(square(sd2F(i)))))*exp(-(square(j-mu2F(i)))/(2*square(sd2F(i))));
              c=Fm(i,j);
 	            if(c2<c) c2=c; 
            }
	 
            for (j=1;j<=12; j++)
              Fm(i,j)=Fm(i,j)/c2;   
            
       }


   Ft = mfexp(log_Ft);

   for (i=1;i<=nyrs; i++)
     {
     for (j=1;j<=12; j++)
          {
           Fm_eff(i,j) = log_Ft(i)*Fm(i,j);    
          }
    }                                      
                   

FUNCTION get_dynamic
 //Recruitments and dynamic of population abundance     
   for (i=1;i<=ncoh; i++)
      {
         id=low(i);  
         N(i,id)=exp(log_Rec(i));
         
	     for (j=low(i)+1;j<=up(i); j++)
         {  
	         yr=year(j-1);
	         edad=edades(i,j-1);
	         mes_cal=meses(j-1);
	         N(i,j)=N(i,j-1)*exp(-1.0*(sel(edad)*Ft(yr)*Fm(yr,mes_cal)+M));
         }
      }
   
 // Generate Predictions
  // Catch in tons
   for (i=1;i<=ncoh; i++)
      for (j=low(i);j<=up(i); j++)
         {
	    yr=year(j);
	    edad=edades(i,j);
	    mes_cal=meses(j);
	    C(i,j)=N(i,j)*((sel(edad)*Ft(yr)*Fm(yr,mes_cal))/(sel(edad)*Ft(yr)*Fm(yr,mes_cal)+M))*(1.0-exp(-1.0*(sel(edad)*Ft(yr)*Fm(yr,mes_cal)+M)));
         } 
  
   for (i=1;i<=nper; i++)
       {
        
        for (j=1;j<=ncoh_alive(i); j++) //number of cohorts alive
          {
	         coh=ncoh_which(i,j);
	         edad=edades(coh,i);
           
            for (k=1;k<=nlen; k++) //
               {
	              CL(i,k)+= C(coh,i)*Age_length(coh,edad,k); //
	              tNatL(i,k)+= N(coh,i)*Age_length(coh,edad,k);
                vNatL(i,k)+= N(coh,i)*sel(edad)*Age_length(coh,edad,k);
                sNatL(i,k)+= N(coh,i)*matur(edad)*Age_length(coh,edad,k);
               }                  
          }
        
        for (k=1;k<=nlen; k++) //multiplication of catch at length by weight at length
          {
            Catch(i)+= CL(i,k)*W_length(i,k);
            CPUE(i)+=vNatL(i,k)*W_length(i,k); // CPUE calculation
            BT(i)+=tNatL(i,k)*W_length(i,k);
            BV(i)+=vNatL(i,k)*W_length(i,k);
            SB(i)+=sNatL(i,k)*W_length(i,k);
            MA(i)+=tNatL(i,k);
          }
          
       for (k=1;k<=9; k++) //multiplication of catch at length by weight at length
          {
            RB(i)+=tNatL(i,k)*W_length(i,k);
          }
      
          pCL(i)=CL(i)/sum(CL(i));      // CL Fishery
   	   }        
   
   CPUE=CPUE*exp(log_q); 

 //*********************************************************************************  
 // CL Reclas
 //*********************************************************************************

   for (l=1;l<=nobs_CL_Reclas; l++) //
       {
  	 i=yrs_CL_Reclas(l); 
         for (j=1;j<=ncoh_alive(i); j++) //number of cohorts
           {
	     coh=ncoh_which(i,j);
	     edad=edades(coh,i);
             for (int k=1;k<=nlen; k++) //
               {
                  vNatL_Reclas_CL(l,k)+= N(coh,i)*sel_Reclas(edad)*Age_length(coh,edad,k);
       	       }                  
           } 
	 CL_Reclas(l)= vNatL_Reclas_CL(l)/sum(vNatL_Reclas_CL(l));    
       }   

 //*********************************************************************************  
 // CL Pelaces
 //*********************************************************************************

  for (l=1;l<=nobs_CL_Pelaces; l++)
       {
  	 i=yrs_CL_Pelaces(l); 
         for (j=1;j<=ncoh_alive(i); j++)
           {
	     coh=ncoh_which(i,j);
	     edad=edades(coh,i);
             for (int k=1;k<=nlen; k++)
               {
                  vNatL_Pelaces_CL(l,k)+= N(coh,i)*sel_Pelaces(edad)*Age_length(coh,edad,k);
       	       }                  
           } 
	 CL_Pelaces(l)= vNatL_Pelaces_CL(l)/sum(vNatL_Pelaces_CL(l));    
       }   
       
 //*********************************************************************************  
 // Index Reclas
 //*********************************************************************************

   for (l=1;l<=nobs_Reclas; l++) //
       {
  	 i=yrs_Reclas(l); 
         for (j=1;j<=ncoh_alive(i); j++) //number of cohorts
           {
	     coh=ncoh_which(i,j);
	     edad=edades(coh,i);
             for (k=1;k<=nlen; k++) //
               {
                  vNatL_Reclas_Index(l,k)+= N(coh,i)*sel_Reclas(edad)*Age_length(coh,edad,k);
       	       }                  
           } 
	
	for (k=1;k<=nlen; k++) //multiplication of catch at length by weight at length
          {
             Index_Reclas(l)+= vNatL_Reclas_Index(l,k)*W_length(i,k);    
          }
       }   
    Index_Reclas=Index_Reclas*exp(Reclas_log_q);
   
 //*********************************************************************************  
 // Index Pelaces
 //*********************************************************************************

   for (l=1;l<=nobs_Pelaces; l++) //
       {
  	 i=yrs_Pelaces(l); 
         for (j=1;j<=ncoh_alive(i); j++) //number of cohorts
           {
	     coh=ncoh_which(i,j);
	     edad=edades(coh,i);
             for (int k=1;k<=nlen; k++) //
               {
                  vNatL_Pelaces_Index(l,k)+= N(coh,i)*sel_Pelaces(edad)*Age_length(coh,edad,k);
       	       }                  
           } 
	
	for (k=1;k<=nlen; k++) //multiplication of catch at length by weight at length
          {
             Index_Pelaces(l)+= vNatL_Pelaces_Index(l,k)*W_length(i,k);    
          }
       }   
    Index_Pelaces=Index_Pelaces*exp(Pelaces_log_q);
 

 
 //*********************************************************************************   
 //Likelihoods
 //*********************************************************************************   

 //likelihood  Catch
  for (i=1;i<=nobs_Catch; i++)
	   {
	     yr=yrs_catch(i);
	     pCatch(i)=Catch(yr); 
	   }
	   LikeType=2;
	   Catch_LL = nlikelihood(obs_Catch, pCatch, cv_Catch, LikeType);
	 
      if(flag_Catch > 0)
         {
	         f+=sum(Catch_LL);           	
   	     }
                                      
    Catch_residuals=(obs_Catch-pCatch)/std_dev(obs_Catch);
  
 // likelihood CPUE
  for (i=1;i<=nobs_CPUE; i++)
	   {
	     yr=yrs_CPUE(i);
	     pCPUE(i)=CPUE(yr); 
	   }
       LikeType=2;
       CPUE_LL = nlikelihood(obs_CPUE,pCPUE,cv_CPUE,LikeType);

        if(flag_CPUE > 0)
         {
	         f+=sum(CPUE_LL);         	
         }
         
   CPUE_residuals=(obs_CPUE-pCPUE)/std_dev(obs_CPUE);      
     
 // Reclas Index Likelihood 
 LikeType=2;
       Reclas_LL = nlikelihood(obs_Reclas,Index_Reclas,cv_Reclas,LikeType);
	 
        if(flag_Reclas > 0)
         {
	          f+=sum(Reclas_LL);         	
   	     }                    
  
  Reclas_residuals=((obs_Reclas-Index_Reclas)/std_dev(obs_Reclas));
  
 // Pelaces Index Likelihood 
 LikeType=2;
     Pelaces_LL=nlikelihood(obs_Pelaces, Index_Pelaces, cv_Pelaces, LikeType);
	 
        if(flag_Pelaces > 0)
         {
	         f+=sum(Pelaces_LL);         	
     	   }

  Pelaces_residuals=(obs_Pelaces-Index_Pelaces)/std_dev(obs_Pelaces);
  
  //*********************************************************************************
  // Catch at length  DATA
  // Fisheries CL
  //*********************************************************************************
  for (i=1;i<=nobs_CL; i++)
	   {
	     yr=yrs_CL(i);
	     predCL(i)=pCL(yr); 
	   }

    int LikeType_CL=12;
     switch (LikeType_CL)
            {
	            case 11:
                 CL_LL = multinomial(obs_CL ,predCL ,nlen , nobs_CL, n_CL);
                   break;
              case 12:
                 CL_LL = nrobust_p(obs_CL ,predCL ,nlen , nobs_CL, n_CL);       
		               break;
	           }

     if(flag_CL > 0)
         {
	         f+=sum(CL_LL);     	
   	     }
   	 //    CL_residuals=obs_CL-predCL;   Parte modificada del codigo
     
     for(i=1;i<=nobs_CL;i++)
     {
       CL_obs_stdev(i)=std_dev(obs_CL(i));
       
       for(j=1;j<=nlen;j++)
       {
       CL_residuals(i,j)=(obs_CL(i,j)-predCL(i,j))/CL_obs_stdev(i);
       }
     }
       
  //*********************************************************************************
  // RECLAS CL
  //*********************************************************************************
  LikeType_CL=12;
     switch (LikeType_CL)
            {
	       case 11:
                 CL_Reclas_LL = multinomial(obs_CL_Reclas,CL_Reclas,nlen,nobs_CL_Reclas,n_CL_Reclas);
                 break;
         case 12:
                 CL_Reclas_LL = nrobust_p(obs_CL_Reclas,CL_Reclas,nlen,nobs_CL_Reclas,n_CL_Reclas);       
		             break;
            }

     if(flag_CL_Reclas > 0)  	// Catch
         {
	         f+=sum(CL_Reclas_LL);         	
     	   }
     	   
     	   CL_Reclas_residuals=obs_CL_Reclas-CL_Reclas;
       
  //*********************************************************************************
  // PELACES CL
  //*********************************************************************************
   LikeType_CL=12;
     switch (LikeType_CL)
            {
	         case 11:
                 CL_Pelaces_LL = multinomial(obs_CL_Pelaces,CL_Pelaces,nlen,nobs_CL_Pelaces,n_CL_Pelaces);
                 break;
                 case 12:
                 CL_Pelaces_LL = nrobust_p(obs_CL_Pelaces,CL_Pelaces,nlen,nobs_CL_Pelaces,n_CL_Pelaces);       
                 break;
	          }
     if(flag_CL_Pelaces > 0)  	// Catch
         {
	         f+=sum(CL_Pelaces_LL);         	
     	   }
       
       
       CL_Pelaces_residuals=obs_CL_Pelaces-CL_Pelaces;
         
  //f=square(dummy);

  f_prof_cl = sum(Catch_LL); 
  
REPORT_SECTION

  report << "Ft \n" << Ft <<endl;
  report << "Fm \n" << Fm <<endl;
  report << "Fm_eff \n" << Fm_eff <<endl;
  report << "Obs_catch \n" << obs_Catch << endl;
  report << "pCatch \n" << pCatch << endl;
  report << "year \n" << year << endl;
  report << "mes \n" << meses << endl;
  report << "obs_Reclas \n" << obs_Reclas << endl;
  report << "Index_Reclas \n" << Index_Reclas << endl;
  report << "obs_Pelaces \n" << obs_Pelaces << endl;
  report << "Index_Pelaces \n" << Index_Pelaces << endl;  
  report << "yrsReclas \n" << yrs_Reclas << endl;
  report << "yrsPelaces \n" << yrs_Pelaces << endl;
  report << "obs_CPUE \n" << obs_CPUE << endl;
  report << "pCPUE \n" << pCPUE << endl;
  report << "yrsCPUE \n" << yrs_CPUE << endl;
  report << "obsCLReclas \n" << obs_CL_Reclas << endl;
  report << "pCLReclas \n" << CL_Reclas << endl;
  report << "yrsCLReclas \n" << yrs_CL_Reclas << endl;
  report << "obsCLPelaces \n" << obs_CL_Pelaces << endl;
  report << "pCLPelaces \n" << CL_Pelaces << endl;
  report << "yrsCLPelaces \n" << yrs_CL_Pelaces << endl;
  report << "selFishery \n" << sel << endl;
  report << "selReclas \n" << sel_Reclas << endl;
  report << "selPelaces \n" << sel_Pelaces << endl;
  report << "maturity \n" << matur << endl;
  report << "edades \n" << edades << endl;
  report << "seasonF \n" << Fm << endl;
  report << "mortality \n" << Ft << endl;



 //=================================================================================  
 // Extra CODES
 //================================================================================= 

FUNCTION get_Likelihoods
   ofstream out1("likelihood.dat",ios::trunc);
   out1<<"**Statistics & Parameters**"<<endl<<endl;
   out1<<*objective_function_value::pobjfun<<endl;
   out1<<objective_function_value::gmax<<endl;   
   out1<<sum(Catch_LL)<<endl;
   out1<<sum(CPUE_LL)<<endl;
   out1<<sum(Reclas_LL)<<endl;
   out1<<sum(Pelaces_LL)<<endl;
   out1<<sum(CL_LL)<<endl;
   out1<<sum(CL_Reclas_LL)<<endl;
   out1<<sum(CL_Pelaces_LL)<<endl;
   out1<<" "<<endl;
   out1<<" "<<endl;
   out1<<Afull<<endl;
   out1<<varL<<endl;
   out1<<Reclas_A_50<<endl;
   out1<<Reclas_ran<<endl;
   out1<<Pelaces_A_50<<endl;
   out1<<Pelaces_ran<<endl;
   out1<<log_q<<endl;
   out1<<Reclas_log_q<<endl;
   out1<<Pelaces_log_q<<endl;
   out1<<mean(log_Rec)<<endl;
   out1<<mean(Ft)<<endl;
   out1<<mean(mu1F)<<endl;
   out1<<mean(mu2F)<<endl;
   out1<<mean(sd1F)<<endl;
   out1<<mean(sd2F)<<endl;
   out1<<mean(pF)<<endl;
   out1.close();
    
FUNCTION get_DFLs
   ofstream out1("DFLs.dat",ios::trunc);
   out1<<predCL<<endl;   
   out1<<obs_CL<<endl;
   out1.close();     

FUNCTION get_residuals
   ofstream out1("residuals_indx.dat",ios::trunc);
   out1<<Catch_residuals<<endl;
   out1<<" "<<endl;
   out1<<CPUE_residuals<<endl;
   out1<<" "<<endl;
   out1<<Reclas_residuals<<endl;
   out1<<" "<<endl;
   out1<<Pelaces_residuals<<endl;
   out1<<" "<<endl;
   out1.close();
   
   ofstream out2("residuals_DFL_catch.dat",ios::trunc);
   out2<<trans(CL_residuals)<<endl;
   /*
   out2<<" "<<endl;
   out2<<trans(CL_Reclas_residuals)<<endl;
   out2<<" "<<endl;
   out2<<trans(CL_Pelaces_residuals)<<endl;
   */
   out2.close();
